﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Cheb
{
    class ChebMethod
    {
        public double[,] u; //искомое решение
        public double[,] b;
        public int n; //число разбиений по икс
        public int m; //число разбиений по игрек
        public double h; // шаг по икс
        public double k; // шаг по игрек
        public int l; // число шагов в цикле метода
        public double[] t; // массив параметров метода

        private double di;

        public ChebMethod()
        {
            n = m = 0;
            h = k = 0;
            setDi();
            u = null;
            b = null;
            l = 0;
            t = null;
        }
        public ChebMethod(int l, int n, double h, int m, double k, double[,] b, double[,] u)
        {
            this.n = n;
            this.m = m;
            this.h = h;
            this.k = k;
            setDi();
            this.u = new double[n + 1, m + 1];
            this.b = new double[(n + 1) , (m + 1)];
            for (int i = 0; i <= n; i++)
            {
                for (int j = 0; j <= m; j++)
                {
                    this.u[i, j] = u[i, j];
                    this.b[i,j] = b[i,j];
                }
            }
            this.l = l;
            setT();

        }
        private void setDi()
        {
            di = -2.0 * (1.0 / (h * h) + 1.0 / (k * k));
        }
        private void setT()
        {
            /*double[,] matr = new double[(n-1) * (m-1), (n-1) * (m-1)];
            double d1 = 1.0 / (h * h);
            double dn = 1.0 / (k * k);

            
            for (int c = 0; c < (n-1)*(m-1) ; c++)
            {
                matr[c, c] = di;
                if (c != 0 && ((c+1)%(n-1) - 1 !=0)) matr[c, c - 1] = d1;
                if (c != (n - 1) * (m - 1) - 1 && ((c + 1) % (n - 1) - 1 != -1)) matr[c, c + 1] = d1;
                if (c > n - 2) matr[c, c - (n - 1)] = dn;
                if (c <= (n - 1) * (m - 1) - n) matr[c, c + (n - 1)] = dn;
            }
            double[] wr =new double[(n-1)*(m-1)];
            double[] wi;
            double[,] wl,wu;
            alglib.rmatrixevd(matr, (n - 1) * (m - 1), 0, out wr, out wi,out wl,out wu);

            int Imax = 0;
            int Imin = 0;
            for (int i = 0; i < (n - 1) * (m - 1); i++)
            {
                if (Math.Abs(wr[Imax]) < Math.Abs(wr[i])) Imax = i;
                if(Math.Abs(wr[Imin]) > Math.Abs(wr[i])) Imin = i;
            }*/

            double lmax = 4.0 * Math.Sin(Math.PI * (n - 1) / (2.0 * n)) * Math.Sin(Math.PI * (n - 1) / (2.0 * n)) / (h * h)
                + 4.0 * Math.Sin(Math.PI * (m - 1) / (2.0 * m)) * Math.Sin(Math.PI * (m - 1) / (2.0 * m)) / (k * k);
            double lmin = 4.0 * Math.Sin(Math.PI / (2.0 * n)) * Math.Sin(Math.PI / (2.0 * n)) / (h * h)
                + 4.0 * Math.Sin(Math.PI / (2.0 * m)) * Math.Sin(Math.PI / (2.0 * m)) / (k * k);
            t = new double[l];
            for (int i = 0; i < l; i++)
            {
                t[i] = -1.0 / ((lmax + lmin) / 2.0 + (lmax - lmin) * Math.Cos(Math.PI * (1 + 2.0 * i) / (2.0 * l)) / 2.0);
                //t[i] = 1.0 / ((Math.Abs(wr[Imax] + wr[Imin])) / 2.0 + (Math.Abs(wr[Imax] - wr[Imin])) * Math.Cos(Math.PI * (1 + 2.0 * i) / (2.0 * l)) / 2.0);
            }

        }
        private double findMaxAbs(double[,] d1, double[,] d2, out int Imax, out int Jmax)
        {
            Imax = 1;
            Jmax = 1;
            for (int i = 1; i < n; i++)
                for (int j = 1; j < m; j++)
                    if(Math.Abs(d1[i,j] - d2[i,j])>Math.Abs(d1[Imax,Jmax] - d2[Imax,Jmax]))
                    {
                        Imax = i;
                        Jmax = j;
                    }
                
            return Math.Abs(d1[Imax,Jmax] - d2[Imax,Jmax]);
        }

        public int method(double e,out double ereal, int count,int Imax,int Jmax)
        {
            
            int s = 0;
            double[,] oldu = new double[n+1,m+1];
            
            double dl = 10;
            
            while (dl > e && s < count)
            {
                s++;                  
                
                for (int i = 0; i < l; i++)
                {
                    for (int p = 0; p <= n; p++)
                        for (int j = 0; j <= m; j++)
                            oldu[p, j] = u[p, j];

                    for (int j = 1; j < m; j++)
                    {
                        for (int p = 1; p < n; p++)
                        {
                            u[p, j] = t[i] * (b[p,j] 
                                - di * oldu[p, j] 
                                - (oldu[p + 1, j] + oldu[p - 1, j]) / (h * h) 
                                - (oldu[p, j + 1] + oldu[p, j - 1]) / (k * k)) + oldu[p, j];
                        }
                    }
                }
                dl = findMaxAbs(u, oldu, out Imax, out Jmax);
            }
            ereal = dl;
            return s;
        }

        




    }
}
